"""
    SatelliteRadiance

Observation operators for satellite radiance observations using radiative transfer models.

This module provides interfaces to:
- CRTM (Community Radiative Transfer Model)
- Simplified radiative transfer for testing
- Forward operator: atmospheric profiles → brightness temperatures
- Adjoint operator: brightness temperature gradients → profile gradients

# Satellite Instruments Supported

## Microwave Sounders
- AMSU-A: 15 channels, temperature sounding
- ATMS: 22 channels, temperature and moisture
- MHS: 5 channels, humidity sounding

## Infrared Sounders
- IASI: 8461 channels, hyperspectral infrared
- AIRS: 2378 channels, hyperspectral infrared
- CrIS: 1305 channels, cross-track infrared

## Imagers
- SEVIRI: 12 channels, geostationary
- ABI: 16 channels, GOES-R series

# Radiative Transfer Physics

The forward model computes brightness temperature Tb from atmospheric profiles:

    Tb(ν) = ∫ τ(ν,p) · B(T(p),ν) · dτ/dp dp + εs · B(Ts,ν) · τs

where:
- ν: wavenumber/frequency
- τ: atmospheric transmittance
- B: Planck function
- T(p): temperature profile
- εs: surface emissivity
- Ts: surface temperature
- τs: surface-to-space transmittance

The adjoint model provides:
- ∂Tb/∂T(p): Jacobian w.r.t. temperature profile
- ∂Tb/∂q(p): Jacobian w.r.t. humidity profile
- ∂Tb/∂εs: Jacobian w.r.t. surface emissivity
"""
module SatelliteRadiance

using LinearAlgebra
using ..Interpolation

export CRTMOperator, SimplifiedRTM
export crtm_forward, crtm_adjoint
export simplified_rtm_forward, simplified_rtm_adjoint
export planck_function, inverse_planck
export compute_transmittance, compute_weighting_function

# Physical constants
const h_planck = 6.62607015e-34  # Planck constant [J·s]
const c_light = 2.99792458e8     # Speed of light [m/s]
const k_boltz = 1.380649e-23     # Boltzmann constant [J/K]
const c1 = 2.0 * h_planck * c_light^2  # First radiation constant
const c2 = h_planck * c_light / k_boltz  # Second radiation constant

"""
    SatelliteInstrument

Enumeration of supported satellite instruments.
"""
@enum SatelliteInstrument begin
    AMSUA
    AMSUB
    ATMS
    MHS
    IASI
    AIRS
    CRIS
    SEVIRI
    ABI
    HIRS
end

"""
    CRTMOperator{T<:AbstractFloat}

Interface to Community Radiative Transfer Model (CRTM).

# Fields
- `instrument::SatelliteInstrument`: Satellite instrument
- `channels::Vector{Int}`: Channel numbers to use
- `sensor_id::String`: CRTM sensor ID
- `coefficient_path::String`: Path to CRTM coefficient files
- `use_crtm::Bool`: Whether CRTM library is available
"""
struct CRTMOperator{T<:AbstractFloat}
    instrument::SatelliteInstrument
    channels::Vector{Int}
    sensor_id::String
    coefficient_path::String
    use_crtm::Bool

    function CRTMOperator{T}(instrument::SatelliteInstrument;
                            channels::Vector{Int}=Int[],
                            sensor_id::String="",
                            coefficient_path::String="",
                            use_crtm::Bool=false) where T
        # Set default sensor_id based on instrument
        sid = isempty(sensor_id) ? get_default_sensor_id(instrument) : sensor_id

        # Set default channels if not provided
        chans = isempty(channels) ? get_default_channels(instrument) : channels

        new{T}(instrument, chans, sid, coefficient_path, use_crtm)
    end
end

CRTMOperator(instrument::SatelliteInstrument; kwargs...) =
    CRTMOperator{Float64}(instrument; kwargs...)

"""
    SimplifiedRTM{T<:AbstractFloat}

Simplified radiative transfer model for testing and development.

Uses simplified physics:
- Absorption-only (no scattering)
- Single gas species (combined H2O/O2/CO2)
- Simplified transmittance calculation
"""
struct SimplifiedRTM{T<:AbstractFloat}
    instrument::SatelliteInstrument
    channels::Vector{Int}
    channel_frequencies::Vector{T}  # GHz for microwave, cm⁻¹ for IR

    function SimplifiedRTM{T}(instrument::SatelliteInstrument;
                             channels::Vector{Int}=Int[]) where T
        chans = isempty(channels) ? get_default_channels(instrument) : channels
        freqs = get_channel_frequencies(instrument, chans)

        new{T}(instrument, chans, freqs)
    end
end

SimplifiedRTM(instrument::SatelliteInstrument; kwargs...) =
    SimplifiedRTM{Float64}(instrument; kwargs...)

# CRTM Interface using Fortran bindings

# Check CRTM availability
const CRTM_AVAILABLE = try
    include("../FortranInterface/crtm_bindings.jl")
    true
catch e
    @warn "CRTM bindings not available, using simplified RTM" exception=e
    false
end

"""
    crtm_forward(operator::CRTMOperator{T},
                t_profile::AbstractVector{T},
                q_profile::AbstractVector{T},
                p_profile::AbstractVector{T};
                surface_temp::T=T(288.0),
                surface_emissivity::T=T(0.95),
                zenith_angle::T=T(0.0),
                azimuth_angle::T=T(0.0)) -> Vector{T}

Forward CRTM operator: atmospheric profiles → brightness temperatures.

# Arguments
- `operator`: CRTM operator configuration
- `t_profile`: Temperature profile [K]
- `q_profile`: Specific humidity profile [kg/kg]
- `p_profile`: Pressure profile [Pa]
- `surface_temp`: Surface temperature [K]
- `surface_emissivity`: Surface emissivity [0-1]
- `zenith_angle`: Satellite zenith angle [degrees]
- `azimuth_angle`: Satellite azimuth angle [degrees]

# Returns
- Vector of brightness temperatures [K] for each channel
"""
function crtm_forward(operator::CRTMOperator{T},
                     t_profile::AbstractVector{T},
                     q_profile::AbstractVector{T},
                     p_profile::AbstractVector{T};
                     surface_temp::T=T(288.0),
                     surface_emissivity::T=T(0.95),
                     zenith_angle::T=T(0.0),
                     azimuth_angle::T=T(0.0)) where T<:AbstractFloat

    nchans = length(operator.channels)
    nlev = length(t_profile)

    if operator.use_crtm && CRTM_AVAILABLE
        # Convert to CRTM format (single profile)
        # CRTM expects (nlev × nprof) matrices
        temp_mat = reshape(Float32.(t_profile), nlev, 1)
        pres_mat = reshape(Float32.(p_profile ./ 100.0), nlev, 1)  # Convert Pa to hPa
        humid_mat = reshape(Float32.(q_profile), nlev, 1)

        tskin = Float32[surface_temp]
        wspd = Float32[5.0]  # Default wind speed
        zenith = Float32[zenith_angle]
        azimuth = Float32[azimuth_angle]

        channels_i32 = Int32.(operator.channels)

        # Call CRTM forward operator
        bt_mat = crtm_forward_julia(temp_mat, pres_mat, humid_mat, tskin, wspd,
                                     zenith, azimuth, channels_i32)

        # Extract and convert to requested type
        bt = Vector{T}(bt_mat[:, 1])
    else
        # Use simplified RTM as fallback
        rtm = SimplifiedRTM{T}(operator.instrument, channels=operator.channels)
        bt = simplified_rtm_forward(rtm, t_profile, q_profile, p_profile,
                                    surface_temp=surface_temp,
                                    surface_emissivity=surface_emissivity)
    end

    return bt
end

"""
    crtm_adjoint(operator::CRTMOperator{T},
                bt_gradient::AbstractVector{T},
                t_profile::AbstractVector{T},
                q_profile::AbstractVector{T},
                p_profile::AbstractVector{T};
                surface_temp::T=T(288.0),
                surface_emissivity::T=T(0.95),
                zenith_angle::T=T(0.0),
                azimuth_angle::T=T(0.0)) -> Tuple{Vector{T}, Vector{T}}

Adjoint CRTM operator: brightness temperature gradients → profile gradients.

# Returns
- `(t_gradient, q_gradient)`: Gradients w.r.t. temperature and humidity profiles
"""
function crtm_adjoint(operator::CRTMOperator{T},
                     bt_gradient::AbstractVector{T},
                     t_profile::AbstractVector{T},
                     q_profile::AbstractVector{T},
                     p_profile::AbstractVector{T};
                     surface_temp::T=T(288.0),
                     surface_emissivity::T=T(0.95),
                     zenith_angle::T=T(0.0),
                     azimuth_angle::T=T(0.0)) where T<:AbstractFloat

    nlev = length(t_profile)
    nchans = length(operator.channels)

    if operator.use_crtm && CRTM_AVAILABLE
        # Convert to CRTM format (single profile)
        temp_mat = reshape(Float32.(t_profile), nlev, 1)
        pres_mat = reshape(Float32.(p_profile ./ 100.0), nlev, 1)  # Convert Pa to hPa
        humid_mat = reshape(Float32.(q_profile), nlev, 1)

        tskin = Float32[surface_temp]
        wspd = Float32[5.0]
        zenith = Float32[zenith_angle]
        azimuth = Float32[azimuth_angle]

        channels_i32 = Int32.(operator.channels)

        # Innovation matrix (nchans × nprof)
        innov_mat = reshape(Float32.(bt_gradient), nchans, 1)

        # Call CRTM adjoint operator
        temp_ad, humid_ad = crtm_adjoint_julia(innov_mat, temp_mat, pres_mat, humid_mat,
                                               tskin, wspd, zenith, azimuth, channels_i32)

        # Extract and convert to requested type
        t_gradient = Vector{T}(temp_ad[:, 1])
        q_gradient = Vector{T}(humid_ad[:, 1])
    else
        # Use simplified RTM adjoint
        rtm = SimplifiedRTM{T}(operator.instrument, channels=operator.channels)
        t_gradient, q_gradient = simplified_rtm_adjoint(rtm, bt_gradient, t_profile, q_profile, p_profile,
                                     surface_temp=surface_temp,
                                     surface_emissivity=surface_emissivity)
    end

    return (t_gradient, q_gradient)
end

# Simplified Radiative Transfer Model

"""
    simplified_rtm_forward(rtm::SimplifiedRTM{T},
                          t_profile::AbstractVector{T},
                          q_profile::AbstractVector{T},
                          p_profile::AbstractVector{T};
                          surface_temp::T=T(288.0),
                          surface_emissivity::T=T(0.95)) -> Vector{T}

Simplified radiative transfer forward model.

Uses a simplified absorption-only model with exponential transmittance.
"""
function simplified_rtm_forward(rtm::SimplifiedRTM{T},
                               t_profile::AbstractVector{T},
                               q_profile::AbstractVector{T},
                               p_profile::AbstractVector{T};
                               surface_temp::T=T(288.0),
                               surface_emissivity::T=T(0.95)) where T<:AbstractFloat

    nlev = length(t_profile)
    nchans = length(rtm.channels)
    bt = zeros(T, nchans)

    for (ic, freq) in enumerate(rtm.channel_frequencies)
        # Compute layer transmittances
        tau = compute_transmittance(freq, t_profile, q_profile, p_profile)

        # Integrate radiative transfer equation
        # Simple upward integration
        radiance = T(0)

        # Surface contribution
        tau_surface = prod(tau)  # Total transmittance to surface
        radiance += surface_emissivity * planck_function(surface_temp, freq) * tau_surface

        # Atmospheric contributions (layer by layer)
        tau_cumulative = T(1)
        for k in 1:nlev
            # Weighting function for this layer
            if k < nlev
                dtau = tau[k] - (k > 1 ? prod(tau[1:k-1]) : T(1))
                layer_emission = planck_function(t_profile[k], freq) * abs(dtau) * tau_cumulative
                radiance += layer_emission
                tau_cumulative *= tau[k]
            end
        end

        # Convert radiance to brightness temperature
        bt[ic] = inverse_planck(radiance, freq)
    end

    return bt
end

"""
    simplified_rtm_adjoint(rtm::SimplifiedRTM{T},
                          bt_gradient::AbstractVector{T},
                          t_profile::AbstractVector{T},
                          q_profile::AbstractVector{T},
                          p_profile::AbstractVector{T};
                          surface_temp::T=T(288.0),
                          surface_emissivity::T=T(0.95)) -> Tuple{Vector{T}, Vector{T}}

Adjoint of simplified radiative transfer model.

Returns gradients with respect to temperature and humidity profiles.
"""
function simplified_rtm_adjoint(rtm::SimplifiedRTM{T},
                               bt_gradient::AbstractVector{T},
                               t_profile::AbstractVector{T},
                               q_profile::AbstractVector{T},
                               p_profile::AbstractVector{T};
                               surface_temp::T=T(288.0),
                               surface_emissivity::T=T(0.95)) where T<:AbstractFloat

    nlev = length(t_profile)
    t_gradient = zeros(T, nlev)
    q_gradient = zeros(T, nlev)

    # Simplified adjoint - distributes gradient to temperature profile
    # Full implementation would compute proper Jacobians
    # This is a placeholder that distributes based on weighting functions

    for (ic, freq) in enumerate(rtm.channel_frequencies)
        grad = bt_gradient[ic]

        # Compute weighting functions
        weights = compute_weighting_function(freq, t_profile, q_profile, p_profile)

        # Distribute gradient according to weighting function
        for k in 1:nlev
            t_gradient[k] += grad * weights[k]
            # Humidity gradient is typically smaller for temperature channels
            q_gradient[k] += grad * weights[k] * T(0.1)
        end
    end

    return (t_gradient, q_gradient)
end

# Utility functions

"""
    planck_function(temperature::T, frequency::T) -> T

Compute Planck radiance for given temperature and frequency.

For microwave (frequency in GHz):
    B(T,ν) = 2kT(ν/c)²

For infrared (frequency in cm⁻¹):
    B(T,ν) = c₁ν³ / (exp(c₂ν/T) - 1)
"""
function planck_function(temperature::T, frequency::T) where T<:AbstractFloat
    # Simplified - assume infrared
    # Convert cm⁻¹ to m⁻¹
    freq_m = frequency * T(100)

    # Planck function
    exponent = c2 * freq_m / temperature
    if exponent > T(80)  # Avoid overflow
        return T(0)
    end

    radiance = c1 * freq_m^3 / (exp(exponent) - T(1))

    return radiance
end

"""
    inverse_planck(radiance::T, frequency::T) -> T

Convert radiance to brightness temperature using inverse Planck function.
"""
function inverse_planck(radiance::T, frequency::T) where T<:AbstractFloat
    if radiance <= T(0)
        return T(0)
    end

    # Simplified - assume infrared
    freq_m = frequency * T(100)

    # Inverse Planck
    ratio = c1 * freq_m^3 / radiance + T(1)
    if ratio <= T(1)
        return T(1000)  # Very hot
    end

    tb = c2 * freq_m / log(ratio)

    return tb
end

"""
    compute_transmittance(frequency::T,
                         t_profile::AbstractVector{T},
                         q_profile::AbstractVector{T},
                         p_profile::AbstractVector{T}) -> Vector{T}

Compute layer transmittances for simplified RTM.

Uses exponential absorption: τ = exp(-k·u)
where k is absorption coefficient and u is absorber amount.
"""
function compute_transmittance(frequency::T,
                              t_profile::AbstractVector{T},
                              q_profile::AbstractVector{T},
                              p_profile::AbstractVector{T}) where T<:AbstractFloat

    nlev = length(t_profile)
    tau = ones(T, nlev)

    # Simplified absorption coefficient (frequency dependent)
    # This is a crude approximation
    k_abs = T(0.001) * frequency  # Increases with frequency

    for k in 1:nlev-1
        # Layer thickness (pressure difference)
        dp = abs(p_profile[k+1] - p_profile[k])

        # Water vapor amount in layer
        q_layer = (q_profile[k] + q_profile[k+1]) / T(2)

        # Optical depth
        optical_depth = k_abs * q_layer * dp / T(1000)

        # Transmittance
        tau[k] = exp(-optical_depth)
    end

    return tau
end

"""
    compute_weighting_function(frequency::T,
                               t_profile::AbstractVector{T},
                               q_profile::AbstractVector{T},
                               p_profile::AbstractVector{T}) -> Vector{T}

Compute channel weighting function for sensitivity to each level.

Weighting function: W(p) = -dτ/dp

Indicates where in the atmosphere the channel is most sensitive.
"""
function compute_weighting_function(frequency::T,
                                   t_profile::AbstractVector{T},
                                   q_profile::AbstractVector{T},
                                   p_profile::AbstractVector{T}) where T<:AbstractFloat

    nlev = length(t_profile)
    weights = zeros(T, nlev)

    # Compute transmittance
    tau = compute_transmittance(frequency, t_profile, q_profile, p_profile)

    # Weighting function is derivative of transmittance
    # Simplified: use finite differences
    for k in 2:nlev-1
        dp = p_profile[k+1] - p_profile[k-1]
        if abs(dp) > eps(T)
            dtau_dp = (tau[k+1] - tau[k-1]) / dp
            weights[k] = -dtau_dp
        end
    end

    # Normalize
    total = sum(abs.(weights))
    if total > eps(T)
        weights ./= total
    end

    return weights
end

# Instrument-specific configuration

"""
    get_default_sensor_id(instrument::SatelliteInstrument) -> String

Get default CRTM sensor ID for instrument.
"""
function get_default_sensor_id(instrument::SatelliteInstrument)
    if instrument == AMSUA
        return "amsua_n19"
    elseif instrument == ATMS
        return "atms_npp"
    elseif instrument == IASI
        return "iasi_metop-a"
    elseif instrument == AIRS
        return "airs_aqua"
    elseif instrument == CRIS
        return "cris_npp"
    else
        return "unknown"
    end
end

"""
    get_default_channels(instrument::SatelliteInstrument) -> Vector{Int}

Get default channel subset for instrument.
"""
function get_default_channels(instrument::SatelliteInstrument)
    if instrument == AMSUA
        return [1, 2, 3, 4, 5, 6, 7, 8]  # Temperature channels
    elseif instrument == ATMS
        return [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    elseif instrument == IASI
        return [16, 38, 49, 51, 55, 57, 59, 61]  # Representative channels
    elseif instrument == AIRS
        return [1, 6, 7, 10, 11, 15, 16, 17]
    else
        return [1, 2, 3, 4, 5]
    end
end

"""
    get_channel_frequencies(instrument::SatelliteInstrument, channels::Vector{Int}) -> Vector{Float64}

Get central frequencies for instrument channels.

Returns GHz for microwave, cm⁻¹ for infrared.
"""
function get_channel_frequencies(instrument::SatelliteInstrument, channels::Vector{Int})
    if instrument == AMSUA
        # AMSU-A frequencies in GHz
        all_freqs = [23.8, 31.4, 50.3, 52.8, 53.596, 54.4, 54.94, 55.5,
                     57.29, 57.29, 57.29, 57.29, 57.29, 57.29, 89.0]
        return all_freqs[channels]

    elseif instrument == IASI
        # IASI: representative wavenumbers in cm⁻¹
        # Full range is 645-2760 cm⁻¹
        representative = [700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0]
        return representative[min.(channels, length(representative))]

    else
        # Default: arbitrary frequencies
        return collect(Float64, 1:length(channels)) .* 100.0
    end
end

end # module SatelliteRadiance
